home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Disc to the Future 2
/
Disc to the Future Part II Programmer's Reference (Wayzata Technology)(6013)(1992).bin
/
MAC
/
THINKC
/
4_0
/
ORBMECHD
/
KEPLER.C
< prev
next >
Wrap
Text File
|
1991-02-18
|
8KB
|
358 lines
/*****************************************************************
Kepler's equation solution algorithms from Battin, Mathematics
and Methods of Astrodynamics, pp192-6.
8 August 1990
*****************************************************************/
#include "orbmech.h"
#include "SANE.h"
#include "math.h"
#define KEPLERREC_SIZE 50
/* globals */
extern DialogPtr KeplerDia;
extern EventRecord gEvent;
extern TEHandle TEH;
extern decform gdecform;
extern int dirty;
struct ParamRec
{
extended value;
Str255 name;
};
extern struct ParamRec gParameter;
struct KeplerRec
{
extended f_1, f_2, time;
extended a, e;
} ;
typedef struct KeplerRec *KeplerPtr;
/*****************************************/
Kepler()
{
Handle inputH, itemHandle;
Str255 f1Str, f2Str, ftStr, aStr, eStr;
int itemType;
Rect itemRect;
Boolean dont_panic = 1, miFlag = 0, auFlag = 0, degFlag = 0;
decimal result;
ControlHandle kmSecBut, miHrBut, radBut, timeBut;
struct KeplerRec *kdata;
Ptr temp;
Str255 q, w, e, r, t, y, u, i;
temp = NewPtr( KEPLERREC_SIZE );
kdata = ( KeplerPtr) temp;
GetDItem( KeplerDia, F_TIME, &itemType, &inputH, &itemRect);
GetIText( inputH, &ftStr );
Pstr_to_extended( &ftStr, &kdata->time, &dont_panic );
GetDItem( KeplerDia, F_1, &itemType, &inputH, &itemRect);
GetIText( inputH, &f1Str );
Pstr_to_extended( &f1Str, &kdata->f_1, &dont_panic );
GetDItem( KeplerDia, F_2, &itemType, &inputH, &itemRect);
GetIText( inputH, &f2Str );
Pstr_to_extended( &f2Str, &kdata->f_2, &dont_panic );
GetDItem( KeplerDia, A_IN, &itemType, &inputH, &itemRect);
GetIText( inputH, &aStr );
Pstr_to_extended( &aStr, &kdata->a, &dont_panic );
GetDItem( KeplerDia, E_IN, &itemType, &inputH, &itemRect);
GetIText( inputH, &eStr );
Pstr_to_extended( &eStr, &kdata->e, &dont_panic );
GetDItem ( KeplerDia, KM_SEC_O, &itemType, &kmSecBut, &itemRect);
GetDItem ( KeplerDia, MI_HR_O, &itemType, &miHrBut, &itemRect);
GetDItem ( KeplerDia, RAD_BUT, &itemType, &radBut, &itemRect);
GetDItem ( KeplerDia, TIME_UNK, &itemType, &timeBut, &itemRect);
if ( dont_panic )
{
if ( ! GetCtlValue ( kmSecBut ) )
if ( GetCtlValue ( miHrBut ) ) {
kdata->a *= MI_CONV;
kdata->time *= HR_CONV;
miFlag = 1;
}
else {
kdata->a *= AU_CONV;
kdata->time *= YR_CONV;
auFlag = 1;
}
if ( ! GetCtlValue ( radBut ) ) {
kdata->f_1 /= DEG_CONV;
kdata->f_2 /= DEG_CONV;
degFlag = 1;
}
if ( GetCtlValue ( timeBut ) ) { /* start time solution */
pStrCopy("\p\rTime between two points ", q);
pStrCopy("\p\r start point: ", w);
pStrCopy("\p\r end point: ", e);
pStrCopy("\p\r a: ", r);
pStrCopy("\p\r e: ", t);
pStrCopy("\p\r primary: ", y);
pStrCopy("\p\r\rTime: ", u);
pStrCopy("\p\r\r", i);
TEInsert( q+1, *q, TEH );
TEInsert( w+1, *w, TEH );
TEInsert( f1Str+1, *f1Str, TEH );
TEInsert( e+1, *e, TEH );
TEInsert( f2Str+1, *f2Str, TEH );
TEInsert( r+1, *r, TEH );
TEInsert( aStr+1, *aStr, TEH );
TEInsert( t+1, *t, TEH );
TEInsert( eStr+1, *eStr, TEH );
TEInsert( y+1, *y, TEH );
TEInsert( gParameter.name+1, *gParameter.name, TEH );
KeplerFindTime( kdata );
if ( miFlag )
kdata->time /= HR_CONV;
if ( auFlag )
kdata->time /= YR_CONV;
if ( kdata->time > 9999.0 || kdata->time < .000001 )
gdecform.style = FLOATDECIMAL;
else
gdecform.style = FIXEDDECIMAL;
num2dec( &gdecform, kdata->time, &result );
dec2str( &gdecform, &result, ftStr );
TEInsert( u+1, *u, TEH );
TEInsert( ftStr+1, *ftStr, TEH );
TEInsert( i+1, *i, TEH );
dirty = 1;
GetDItem( KeplerDia, F_TIME, &itemType, &itemHandle, &itemRect);
SetIText( itemHandle, ftStr );
ShowSelect();
}
else { /**** find position block ***/
pStrCopy("\p\rPosition after elapsed time", q);
pStrCopy("\p\r start point: ", w);
pStrCopy("\p\r time: ", e);
pStrCopy("\p\r a: ", r);
pStrCopy("\p\r e: ", t);
pStrCopy("\p\r primary: ", y);
pStrCopy("\p\r\rFinal position: ", u);
pStrCopy("\p\r\r", i);
TEInsert( q+1, *q, TEH );
TEInsert( w+1, *w, TEH );
TEInsert( f1Str+1, *f1Str, TEH );
TEInsert( e+1, *e, TEH );
TEInsert( ftStr+1, *ftStr, TEH );
TEInsert( r+1, *r, TEH );
TEInsert( aStr+1, *aStr, TEH );
TEInsert( t+1, *t, TEH );
TEInsert( eStr+1, *eStr, TEH );
TEInsert( y+1, *y, TEH );
TEInsert( gParameter.name+1, *gParameter.name, TEH );
KeplerFindPosition( kdata );
if ( degFlag )
kdata->f_2 *= DEG_CONV;
if ( kdata->f_2 > 9999.0 || kdata->time < .000001 )
gdecform.style = FLOATDECIMAL;
else
gdecform.style = FIXEDDECIMAL;
num2dec( &gdecform, kdata->f_2, &result );
dec2str( &gdecform, &result, f2Str );
TEInsert( u+1, *u, TEH );
TEInsert( f2Str+1, *f2Str, TEH );
TEInsert( i+1, *i, TEH );
dirty = 1;
GetDItem( KeplerDia, F_2, &itemType, &itemHandle, &itemRect);
SetIText( itemHandle, f2Str );
ShowSelect();
}
}
}
/*****************************/
extended trueToEccentric( f, e )
extended *f, *e;
{
extended ecc;
ecc = 2 * atan( sqrt( (1 - *e) / (1 + *e) ) * tan( 0.5 * *f ) );
return (ecc);
}
/*****************************/
extended eccentricToTrue( ecc, e )
extended *ecc, *e;
{
extended f;
f = 2 * atan( sqrt( (1 + *e) / (1 - *e) ) * tan( 0.5 * *ecc ) );
return (f);
}
/*****************************/
KeplerFindTime( kdata )
struct KeplerRec *kdata;
{
extended f1, f2, k, a, period, ecc1, ecc2, m1, m2;
extended t1, t2, e, temptime, tf1, tf2;
f1 = kdata->f_1;
f2 = kdata->f_2;
tf1 = kdata->f_1;
tf2 = kdata->f_2;
a = kdata->a;
e = kdata->e;
k = sqrt( a * a * a / gParameter.value );
period = 2 * PI_VALUE * k;
if ( tf1 > ( PI_VALUE - 0.00001 ) && tf1 < ( PI_VALUE + 0.00001 ) )
t1 = period / 2.0;
else {
if ( tf1 >= ( PI_VALUE - 0.00001 ) )
tf1 = 2 * PI_VALUE - tf1;
ecc1 = trueToEccentric( &tf1, &e );
t1 = k * ( ecc1 - e * sin( ecc1 ) );
}
if ( tf2 > ( PI_VALUE - 0.00001 ) && tf2 < ( PI_VALUE + 0.00001 ) )
t2 = period / 2.0;
else {
if ( tf2 >= ( PI_VALUE - 0.00001 ) )
tf2 = 2 * PI_VALUE - tf2;
ecc2 = trueToEccentric( &tf2, &e );
t2 = k * ( ecc2 - e * sin( ecc2 ) );
}
temptime = 0;
if ( ( f1 > f2 ) || ( f2 > PI_VALUE && f1 < PI_VALUE ) )
temptime += period;
if ( f1 < PI_VALUE )
t1 *= -1.0;
if ( f2 > PI_VALUE )
t2 *= -1.0;
temptime += t1 + t2;
kdata->time = temptime;
}
/***************************************/
KeplerFindPosition( kdata )
struct KeplerRec *kdata;
{
extended k, a, period, ecc1, m1, t1, ttot, time, f1;
extended f2, ecc2, m2, t2, old_ecc, e;
int converged;
f1 = kdata->f_1;
time = kdata->time;
a = kdata->a;
e = kdata->e;
k = sqrt( a * a * a / gParameter.value );
period = 2 * PI_VALUE * k;
if ( f1 > ( PI_VALUE - 0.00001 ) && f1 < ( PI_VALUE + 0.00001 ) )
t1 = period / 2.0;
else {
if ( f1 >= ( PI_VALUE - 0.00001 ) )
f1 = 2 * PI_VALUE - f1;
ecc1 = trueToEccentric( &f1, &e );
t1 = k * ( ecc1 - e * sin( ecc1 ) );
}
ttot = t1 + time;
while ( ttot > period )
ttot -= period;
converged = 0;
old_ecc = 0;
ecc2 = 0;
m2 = ttot / k;
while ( ! converged ) {
ecc2 = m2 + e * sin( ecc2 );
converged = checkConvergence( ecc2, old_ecc );
old_ecc = ecc2;
}
f2 = eccentricToTrue( &ecc2, &e );
if ( f2 < 0.0 )
f2 += 2 * PI_VALUE;
kdata->f_2 = f2;
}